home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / hyperg_0F1.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  5.1 KB  |  179 lines
  1. /* specfunc/hyperg_0F1.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_gamma.h>
  27. #include <gsl/gsl_sf_bessel.h>
  28. #include <gsl/gsl_sf_hyperg.h>
  29.  
  30. #include "error.h"
  31.  
  32. #define locEPS  (1000.0*GSL_DBL_EPSILON)
  33.  
  34.  
  35. /* Evaluate bessel_I(nu, x), allowing nu < 0.
  36.  * This is fine here because we do not not allow
  37.  * nu to be a negative integer.
  38.  * x > 0.
  39.  */
  40. static
  41. int
  42. hyperg_0F1_bessel_I(const double nu, const double x, gsl_sf_result * result)
  43. {
  44.   if(x > GSL_LOG_DBL_MAX) {
  45.     OVERFLOW_ERROR(result);
  46.   }
  47.  
  48.   if(nu < 0.0) { 
  49.     const double anu = -nu;
  50.     const double s   = 2.0/M_PI * sin(anu*M_PI);
  51.     const double ex  = exp(x);
  52.     gsl_sf_result I;
  53.     gsl_sf_result K;
  54.     int stat_I = gsl_sf_bessel_Inu_scaled_e(anu, x, &I);
  55.     int stat_K = gsl_sf_bessel_Knu_scaled_e(anu, x, &K);
  56.     result->val  = ex * I.val + s * (K.val / ex);
  57.     result->err  = ex * I.err + fabs(s * K.err/ex);
  58.     result->err += fabs(s * (K.val/ex)) * GSL_DBL_EPSILON * anu * M_PI;
  59.     return GSL_ERROR_SELECT_2(stat_K, stat_I);
  60.   }
  61.   else {
  62.     const double ex  = exp(x);
  63.     gsl_sf_result I;
  64.     int stat_I = gsl_sf_bessel_Inu_scaled_e(nu, x, &I);
  65.     result->val = ex * I.val;
  66.     result->err = ex * I.err + GSL_DBL_EPSILON * fabs(result->val);
  67.     return stat_I;
  68.   }
  69. }
  70.  
  71.  
  72. /* Evaluate bessel_J(nu, x), allowing nu < 0.
  73.  * This is fine here because we do not not allow
  74.  * nu to be a negative integer.
  75.  * x > 0.
  76.  */
  77. static
  78. int
  79. hyperg_0F1_bessel_J(const double nu, const double x, gsl_sf_result * result)
  80. {
  81.   if(nu < 0.0) { 
  82.     const double anu = -nu;
  83.     const double s   = sin(anu*M_PI);
  84.     const double c   = cos(anu*M_PI);
  85.     gsl_sf_result J;
  86.     gsl_sf_result Y;
  87.     int stat_J = gsl_sf_bessel_Jnu_e(anu, x, &J);
  88.     int stat_Y = gsl_sf_bessel_Ynu_e(anu, x, &Y);
  89.     result->val  = c * J.val - s * Y.val;
  90.     result->err  = fabs(c * J.err) + fabs(s * Y.err);
  91.     result->err += fabs(anu * M_PI) * GSL_DBL_EPSILON * fabs(J.val + Y.val);
  92.     return GSL_ERROR_SELECT_2(stat_Y, stat_J);
  93.   }
  94.   else {
  95.     return gsl_sf_bessel_Jnu_e(nu, x, result);
  96.   }
  97. }
  98.  
  99.  
  100. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  101.  
  102. int
  103. gsl_sf_hyperg_0F1_e(double c, double x, gsl_sf_result * result)
  104. {
  105.   const double rintc = floor(c + 0.5);
  106.   const int c_neg_integer = (c < 0.0 && fabs(c - rintc) < locEPS);
  107.  
  108.   /* CHECK_POINTER(result) */
  109.  
  110.   if(c == 0.0 || c_neg_integer) {
  111.     DOMAIN_ERROR(result);
  112.   }
  113.   else if(x < 0.0) {
  114.     gsl_sf_result Jcm1;
  115.     gsl_sf_result lg_c;
  116.     double sgn;
  117.     int stat_g = gsl_sf_lngamma_sgn_e(c, &lg_c, &sgn);
  118.     int stat_J = hyperg_0F1_bessel_J(c-1.0, 2.0*sqrt(-x), &Jcm1);
  119.     if(stat_g != GSL_SUCCESS) {
  120.       result->val = 0.0;
  121.       result->err = 0.0;
  122.       return stat_g;
  123.     }
  124.     else if(Jcm1.val == 0.0) {
  125.       result->val = 0.0;
  126.       result->err = 0.0;
  127.       return stat_J;
  128.     }
  129.     else {
  130.       const double tl = log(-x)*0.5*(1.0-c);
  131.       double ln_pre_val = lg_c.val + tl;
  132.       double ln_pre_err = lg_c.err + 2.0 * GSL_DBL_EPSILON * fabs(tl);
  133.       return gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
  134.                                       sgn*Jcm1.val, Jcm1.err,
  135.                       result);
  136.     }
  137.   }
  138.   else if(x == 0.0) {
  139.     result->val = 1.0;
  140.     result->err = 1.0;
  141.     return GSL_SUCCESS;
  142.   }
  143.   else {
  144.     gsl_sf_result Icm1;
  145.     gsl_sf_result lg_c;
  146.     double sgn;
  147.     int stat_g = gsl_sf_lngamma_sgn_e(c, &lg_c, &sgn);
  148.     int stat_I = hyperg_0F1_bessel_I(c-1.0, 2.0*sqrt(x), &Icm1);
  149.     if(stat_g != GSL_SUCCESS) {
  150.       result->val = 0.0;
  151.       result->err = 0.0;
  152.       return stat_g;
  153.     }
  154.     else if(Icm1.val == 0.0) {
  155.       result->val = 0.0;
  156.       result->err = 0.0;
  157.       return stat_I;
  158.     }
  159.     else {
  160.       const double tl = log(x)*0.5*(1.0-c);
  161.       const double ln_pre_val = lg_c.val + tl;
  162.       const double ln_pre_err = lg_c.err + 2.0 * GSL_DBL_EPSILON * fabs(tl);
  163.       return gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
  164.                                       sgn*Icm1.val, Icm1.err,
  165.                       result);
  166.     }
  167.   }
  168. }
  169.  
  170.  
  171. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  172.  
  173. #include "eval.h"
  174.  
  175. double gsl_sf_hyperg_0F1(const double c, const double x)
  176. {
  177.   EVAL_RESULT(gsl_sf_hyperg_0F1_e(c, x, &result));
  178. }
  179.